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Abstract 

Most theories of evolutionary diversification are based on equilibrium assumptions: they are either based 
on optimality arguments involving static fitness landscapes, or they assume that populations first evolve to 
an equilibrium state before diversification occurs, as exemplified by fhe concepf of evolufionary branching 
poinfs in adaptive dynamics fheory. Recenf resulfs indicafe fhaf adapfive dynamics may offen nof con¬ 
verge fo equilibrium poinfs and insfead generafe complicafed frajecfories if evolution lakes place in high¬ 
dimensional phenofype spaces. Even fhough some analytical resulfs on diversificalion in complex pheno¬ 
type spaces are available, fo sfudy fhis problem in general we need fo reconsfrucf individual-based models 
from fhe adapfive dynamics generaling fhe non-equilibrium dynamics. Here we firsl provide a mefhod fo 
conslrucl individual-based models such fhaf fhey faifhfully reproduce fhe given adapfive dynamics aflraclor 
wilhoul diversificalion. We Ihen show fhaf a propensily fo diversify can by inlroduced by adding Gaussian 
compefilion terms that generate frequency dependence while still preserving the same adaptive dynamics. 
For sufficiently strong competition, the disruptive selection generated by frequency-dependence overcomes 
the directional evolution along the selection gradient and leads to diversification in phenotypic directions 
that are orthogonal to the selection gradient. 
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I. INTRODUCTION 


Understanding the origin of biologieal diversity is one of the most fundamental problems in 
evolutionary biology. Traditional theories of diversifieation are based on statie fitness landseapes 
and geographie isolation [|5lll5l. However, it has been realized that eeologieal interaetions lead¬ 
ing to frequeney-dependent seleetion and dynamieally ehanging fitness landseapes ean generate 
adaptive diversifieation, a proeess that oeeurs without geographie isolation and requires eeologieal 
eontaet between the newly emerging speeies iTTi rTOlfT^ . 

There is a substantial amount of empirieal evidenee for adaptive diversifieation in sympatry 
ni 12 HH [191 l22l l26] - [3T|| , and a plethora of different models have shown the theoretieal feasibil¬ 
ity of adaptive diversifieation ( ifTOl . see also the eomprehensive list of papers on Eva Kisdi’s site 
mathstat.helsinki.fi/~kisdi/). In partieular, the framework of adaptive dynamies has been exten¬ 
sively used in this eontext, beeause it is ideally suited to deseribe evolution on fitness landseapes 
that are ehanging dynamieally due to frequeney-dependent interaetions, and to identify eonditions 
that are eondueive to adaptive diversifieation [[TOl [I6l |25l . In essenee, adaptive dynamies is a for¬ 
malism for deriving evolutionary trajeetories in phenotype spaee. In most oases, this phenotype 
spaoe is one-dimensional and represents soalar traits suoh as body size or other morphologioal or 
behavioural features. In this ease, the eonditions for adaptive diversifieation are very well under¬ 
stood and are based on the oonoept of evolutionary branohing points. These are points in pheno¬ 
type spaoe with two essential properties: first, they are attraetors for the evolutionary dynamies, 
and seeond, they are fitness minima. This may seem oontradiotory at first but beeomes intuitively 
appealing onoe one takes into aeeount that fitness landseapes are dynamio. Thus, as long as the 
evolutionary trajeotory is away from the branehing point, the fitness landsoape determining the 
trajeotory has a minimum so that the ourrent phenotypio position of the evolving population is on 
one side of this minimum. As the population elimbs away from the fitness minimum, the fitness 
landsoape ohanges in suoh a way that the minimum eventually oatehes up with the trajeotory, at 
whieh point the dynamies has reaehed its equilibrium, while the population sits on a fitness mini¬ 
mum. Beeause of this, phenotypie mutants on either side of the eurrent position of the population 
ean invade, resulting in adaptive diversifieation. For example, if this soenario is modelled with 
individual-based models, in whieh a population is represented as a oloud of points in phenotype 
spaee, then this oloud first oonverges towards the position in phenotype spaee oorresponding to the 
evolutionary branehing point, and subsequently splits into two separate elouds that diverge from 
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each other, thus yielding and elementary representation of evolutionary diversification. The exis¬ 
tence of evolutionary branching points is a robust feature of adaptive dynamics in low-dimensional 
phenotype spaces and has been demonstrated in many different settings IITOll . 

In the majority of models studied to date, adaptive diversification occurs according to the 
scheme just described: first the evolutionary trajectory converges to a branching point, and then 
the population splits into diverging and coexisting phenotypic branches. For example, it has been 
shown that, conditional on convergence to an evolutionary equilibrium, the likelihood of evolu¬ 
tionary branching increases with the dimension of phenotype space [lS|T3l[32l. However, it has 
recently also been shown that in high-dimensional phenotype spaces, the adaptive dynamics may 
never converge to an equilibrium in the first place ifT^ . In fact, in a general class of competi¬ 
tion models, the likelihood that adaptive dynamics are chaotic approaches 1 if the dimension of 
phenotype space is high enough ifT^ . In particular, such dynamics would never converge to an 
evolutionary branching point. A resolution of this apparent paradox has been found in [[2T]|, which 
shows analytically that sufficiently strong disruptive selection can “overcome” a small but non¬ 
zero selection gradient and produce sustained diversification in the direction perpendicular to the 
selection gradient. The analysis in [l^ is local in the sense that it describes the initial phase of 
evolutionary branching when the resident is moving along an assumed and constant gradient in 
one phenotypic direction. Here we study the question of diversification in non-stationary resident 
populations when the selection gradient along which the resident is evolving is not assumed, but 
is a result of complex interactions in high-dimensional phenotype space. In particular, the trajec¬ 
tory of the resident is not simply given by a constant gradient, but can exhibit more complicated 
dynamics, such as cycles and chaos. 

In the following, we first briefly recall how the conditions for diversification of [l2T1l apply in 
scenarios with complicated evolutionary dynamics of the resident. To study the full dynamics of 
diversification in such scenarios, we then reconstruct individual-based models from an adaptive 
dynamics model that generates complicated trajectories. It is known [HI IH [8l that adaptive dy¬ 
namics is a particular large-population limit of underlying stochastic, individual-based processes. 
What seems to be less appreciated is that there are in fact very many different individual-based 
models that yield the same adaptive dynamics. However, these different individual-based mod¬ 
els differ in their propensity to diversify. Here we show that given an adaptive dynamics model, 
there is a particular “minimal” individual-based model that never diversifies and that produces 
qualitatively the same trajectories as the given adaptive dynamics model. We then show how. 
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in light of the conditions for diversification in [l2T1l . this minimal model can be altered to gener¬ 
ate adaptive diversification, yet still have the same adaptive dynamics limit. Similar results hold 
for partial differential equation (PDE) models, which represent an intermediate limit between the 
individual-based and the adaptive dynamics [[31IH. Our work appears to be one of the first attempts 
to study non-equilibrium evolutionary dynamics of complex, high-dimensional phenotypes using 
individual-based and PDE models. 


II. ADAPTIVE DYNAMICS AND DIVERSIFICATION 


As in [|T2l, we study a general class of models for frequency-dependent competition in which 
ecological interactions are defined by d-dimensional phenotypes, where d > 1. Eor example, one 
can imagine that a high-dimensional phenotype of individuals is given by the efficiencies of several 
metabolic pathways, or my multiple morphological characteristics. The ecological interactions are 
determined by a competition kernel a(x, y) and by a carrying capacity K (x), where x, y G are 
the d-dimensional phenotypes of competing individuals. The competition kernel a measures the 
competitive impact that an individual of phenotype x has on an individual of phenotype y, and in 
the sequel we always assume that a(x, x) = 1 for all x. Assuming logistic ecological dynamics, 
K (x) is then the equilibrium density of a population that is monomorphic for phenotype x. The 
per capita growth rate of a rare mutant phenotype y in the resident x is then given by 


/(x,y) 


a(x,y)A:(x) 

K{y) 


( 1 ) 


(see [TTOlfT^ for more details). The function /(x, y) is the invasion fitness from which the adaptive 
dynamics is derived based on the selection gradients 
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The selection gradients define a system of differential equations on phenotype space which is 
given by 


dx 

dt 


M • s(x). 


(3) 


Here M is the mutational variance-covariance matrix, which generally influences the speed and 
direction of evolution. For simplicity, we assume here that this matrix is the unit matrix. For more 
details on the derivation of the adaptive dynamics Q we refer to a large body of primary literature 
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(e.g. lISUTOl [m |23l). We note that the adaptive dynamies Q ean be derived analytieally as a 
large-population limit of an underlying stoehastie, individual-based model that is again defined by 
the eompetition kernel a(x, y) and the earrying eapaeity -ff (x) [|3l|4l[8l. 

The system of ODEs Q describes the trajectory of an evolving monomorphic population in 
phenotype space Typically, the goal of an adaptive dynamics analysis is to find equilibrium 
attractors for these trajectories, which are called singular points. These are points x* in phenotype 
space at which the selection gradient vanishes, s(x*) = 0. However, in general the trajectories may 
also exhibit more complicated dynamics, such as limit cycles or chaos. In fact, [fl^ have argued 
that in high-dimensional phenotype spaces, almost all systems of the form Q will be chaotic. 

For scalar traits, i.e., when the dimension of phenotype space is one, the adaptive dynamics 
typically converges to a singular point x*, and it is well known that these points are often evo¬ 
lutionary branching points, i.e., starting points for adaptive diversification. This can be seen by 
expanding the invasion fitness function with respect to the mutant coordinate y to second order. 


f{x,y) =f{x,x) + 


df{x,y) 


dy 


{y-x) + 


d^f{x,y) 


y=x 


dy^ 


{y - xY 


(4) 


y=x 


The first term on the right hand side, f{x,x), is zero for all x since the resident population is 
assumed to be at the ecological equilibrium. The coefficient of the linear term is the selection 
gradient, which vanishes at a singular point by definition. In the neighbourhood of singular points, 
df{x, y)/dy\y=x —)■ 0, and the second order term in eq. Q becomes significant. In particular, if 
the singular point is a maximum of the invasion fitness, i. e. d'^f{x*,y)/dy‘^\y=x* < 0, no nearby 
mutants can invade and the population remains monomorphic. However, when the invasion fitness 
has a minimum at the singular point, i. e. when d‘^f{x*, y)/ dyYy=x* > 0, the stationary point is 
an evolutionary branching point. Evolutionary branching has been studied extensively for scalar 
traits ([[IOl[IIl[l6l|2l). 

An extension of these arguments to high-dimensional phenotype spaces reveals two interesting 
trends. On the one hand, the conditions for evolutionary branching at a singular point, i.e., at an 
equilibrium of the adaptive dynamics, are progressively easier to satisfy for increasing dimension¬ 
ality [|^[T3l|32l. On the other hand, and as mentioned before, the fraction of trajectories that con¬ 
verge to a singular point decreases and essentially vanishes for very high-dimensional phenotypes, 
as more and more trajectories become chaotic iTT^ . Thus, adaptive diversification by means of 
convergence to evolutionary branching points and subsequent diversification becomes less likely 
in high-dimensional phenotype spaces, and the possibility of non-equilibrium adaptive dynamics 
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makes it necessary to study diversification not just from evolutionary branching points, but more 
generally from any point on an evolutionary trajectory. Indeed, [|^ have studied this problem 
analytically under the assumption that the resident is under directional selection in a particular 
phenotypic direction. 

To see how the analysis of [l2T]l applies in cases with complicated evolutionary trajectories, 
consider the multidimensional generalization of the expansion of the invasion fitness Q, 


/(x,y) =/(x,x) + 


a/(x,y) 9/(x,y) 


dyi 


dyd 


y=x 


(y-x) + -(y-xfH(x)(y-x), (5) 


where H(x) is the Hessian matrix of second derivatives of the invasion fitness function with respect 
to the mutant trait value y, evaluated at the resident trait value x: 




Hlxl = 
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( 6 ) 


^=x/ 


dvddyd^^^'’ y^ 

As in the one-dimensional case, the selection gradient, which determines the evolutionary tra¬ 
jectory, is the linear term in Q. However, the linear term is not relevant for the remaining d — 1 
dimensions that are orthogonal to the selection gradient. In this orthogonal subspace, the quadratic 
terms determine whether invasion is possible. Specifically, the interplay between the eigenvalues 
of the Hessian matrix H(x) and the selection gradient determine the curvature of the invasion fit¬ 
ness in the directions that are orthogonal to the direction in which the population is evolving [l^ . 
In particular, if H(x) is negative definite, then the fitness function, viewed as a function on the 
subspace that is orthogonal to the selection gradient, has a maximum at the current resident, and 
only mutant phenotypes that have a component in the direction of the selection gradient can in¬ 
vade. Thus, in this case no diversification is expected, and instead the evolving population simply 
follows the trajectory determined by the selection gradient. However, if the Hessian matrix has an 
eigenvector that has a positive eigenvalue and a non-zero projection onto the orthogonal subspace, 
then there are directions in phenotype space that are orthogonal to the selection gradient and along 
which the diversification is possible. The precise conditions for diversification in terms of the 
curvature of the invasion fitness in the direction orthogonal to the selection gradient and the mag¬ 
nitude of the selection gradient have been investigated in [|211. Roughly speaking, diversification 
occurs if the eigenvalues of the orthogonal Hessian are positive enough. 

It is important to note that if the population evolves on a non-equilibrium trajectory due to the 
first order term in Q, then the direction of the selection gradient continually changes, and hence 


6 














so does the subspaee that is orthogonal to the seleetion gradient. In partieular, phenotypie diree- 
tions along whieh the invasion fitness has a minimum may ehange over time, whieh may impede 
diversifieation, as the eonditions for diversifieation may be satisfied only temporary. Neverthe¬ 
less, it seems elear that in prineiple, the rather restrietive eonditions for adaptive diversifieation in 
sealar traits, requiring the adaptive dynamies to eonverge to a branehing point, ean easily be eir- 
eumvented in high-dimensional phenotype spaees, in whieh diversifieation ean oeeur in direetions 
that are orthogonal to the eurrent direetion of evolution. Below we use individual-based models 
to eonfirm that the essential positivity of the loeal eurvature of the Hessian of the invasion fitness 
is indeed neeessary and suffieient for adaptive diversifieation during non-equilibrium evolutionary 
dynamies. 


III. RECONSTRUCTING INDIVIDUAL-BASED MODELS WITH DIEEERENT PATTERNS OE 
DIVERSIEICATION 

Our main goal is to eonstruet individual-based models that eorrespond to a given adaptive 
dynamies Q. In individual-based models, a population is represented by a “eloud” of points 
that moves through phenotype spaee over time. The points in the eloud are the individuals, eaeh 
of whieh is assigned a birth rate and a death rate at any given point in time. The evolutionary 
dynamies unfolds through an algorithm aeeording to whieh individuals with higher birth rates 
are more likely to give birth, and individuals with higher death rates are more likely to die. For 
birth events, the simplest assumption is asexual reproduetion, and when individuals give birth, the 
offspring has a phenotype that is similar to the phenotype of the parent. There are various ways of 
implementing sueh algorithms, and the one used below is deseribed in the Supplementary Online 
Materials (SOM). We refer to the literature for a more detailed deseription of sueh algorithms (e.g. 
E |4l [El)- In sueh individual-based models, a monomorphie population eorresponds to a single 
eompaet eloud, and diversifieation oeeurs when an evolving population splits into two elusters 
in phenotype spaee. Thus, after diversifieation, there are multiple eoexisting elouds of points in 
phenotype spaee, eaeh representing a deseendent lineage. 

An alternative way to deseribe the evolving population is by using partial differential equations 
(PDF’s) for the dynamies of population distributions in phenotype spaee (see SOM). In that ease, 
diversifieation eorresponds to the formation of multiple modes in the evolving distribution, in faet, 
these deterministie models are themselves large-populations limits of the individual-based models. 
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but they are more general than the adaptive dynamies limit in the sense that fewer assumptions are 
neeessary to derive the PDE limit. This has been nieely deseribed in [Q. Due to eomputational 
limitations of numerieal solutions of PDE’s, solving the PDE models is eurrently only feasible in 
relatively low-dimensional phenotype spaees. 

When eonstrueting either-individual based models or PDE models that eorrespond to a given 
adaptive dynamies Q, a fundamental problem arises. The adaptive dynamies Q is determined 
by the seleetion gradient Q, whieh is in turn determined by the derivatives of the invasion fitness 
funetion Q, and henee by the derivatives of the eompetition kernel and the earrying eapaeity. Be- 
eause taking derivatives does obviously not yield a one-to-one eorrespondenee, there are infinitely 
many invasion fitness funetions with different eompetition kernels and earrying eapaeities that 
have the same derivatives, and henee yield the same adaptive dynamies. Sinee both individual- 
based models and PDE models require knowledge of the full eompetition kernel and earrying ea¬ 
paeity (see SOM), this implies that there are infinitely many different sueh models that eorrespond 
to the same adaptive dynamies Q. Importantly, while yielding the same adaptive dynamies, the 
different eompetition kernels and earrying eapaeities produee different patterns of diversifieation. 


A. The “minimal” individual-based model 


To illustrate the above points, we reeonstruet individual-based models from a given adaptive 
dynamies model by writing the adaptive dynamies equations (|^|^ as 

dxi c)a(x, y) 

|y=x 


dt 


dVi 


c)ln[/f(x)] 

H -^= U!i{x) + Ui{x), 


(7) 


where we denote by tej(x) and Mi(x) the terms eoming from the eompetition kernel and the earry¬ 
ing eapaeity, respeetively. 

We first note that the part of the seleetion gradient that is due to the earrying eapaeity funetion 
(i.e., the Ui (x) terms) essentially generate a hill elimbing proeess on the funetion K (x). As in lfT2ll 
we assume that this hill-elimbing proeess is of the simplest possible form by assuming that the 
earrying eapaeity funetion is given by the unimodal funetion 


K (x) = Ko exp 



( 8 ) 


This implies that the terms Mi(x) = —xf eorrespond to a stabilizing eomponent of seleetion, with 
X = 0 being the optimal phenotype, beeause this phenotype eorresponds to the highest earrying 
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capacity. 

Under these assumptions, we first show that regardless of the speeifie form of the Wj(x)-terms 
in (|7]), there is a “minimal” eompetition kernel for which the Hessian matrix of the invasion fitness 
is always negative definite and henee the evolving population should remain monomorphic in the 
individual-based model. 

Specifically, we define the minimal eompetition kernel as 

d 

a(x,y) = exp (9) 

i=l 

where the Wi(x) are taken from the adaptive dynamics Q. It is easy to see that the adaptive 
dynamies corresponding to this eompetition kernel, as well as to the earrying capaeity is given 
by Q as desired. Moreover, with this eompetition kernel, we ealeulate the elements of the Hessian 
(7) of the invasion fitness function (2) as 


H(x)jj = - M;i(x)tej(x) + Wi(x)uj(x) + Wj(x)ui(x) - Ui(x)uj(x) + 


1 d^K(x) 
K(x) dxidxj 


( 10 ) 


This matrix and can be written more sueeinetly as 


H(x) = - 


^ Wi — Ui^ 


\Wd — UdJ 


Wi - Ui, ..., Wd - Ud 


Hxfx) 


K{x) ’ 


( 11 ) 


where Hk(x) is the Hessian of the earrying eapaeity, whieh is by assumption negative definite 
(and where we have omitted the argument x from the terms Wi{x) and Ui{x)). 

The first term in ( [TT] ) is a rank 1 matrix and its only non-zero eigenvalue A is negative: 

d 

A =-^[tei(x)-Mj(x)]^ (12) 


As eonsequenee, with the minimal eompetition kernel the Hessian of the invasion fitness fune- 
tion minimal eompetition kernel is always negative definite, independent of the eurrent resident 
phenotype, and henee diversification is not expected anywhere along the trajeetory of the adaptive 
dynamies. 

As an example, consider the adaptive dynamies studied in f[T^ . 

d d 

bijXj -f ttijkXjXk - x\, i = 1,..., d. (13) 


dxi 

dt 


i=i 


j,k=l 
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where bij and aijk are arbitrary eoeffieients. The eorresponding minimal eompetition kernel is 
derived from these eoeffieients as 


d 


d 



(14) 


while the earrying eapaeity is given by ([^ as before. The predietion is that no matter what adaptive 
dynamies the eoeffieients bij and aijk generate, the eorresponding individual-based model obtained 
from the eompetition kernel ( [T4| ) and the earrying eapaeity ([^ do not show diversifieation and 
produee elouds of points that move roughly along the same trajeetories as the adaptive dynamies 
model ( fT3] ). 

Our extensive simulations of individual-based models for many different ehoiees of the eoef¬ 
fieients bij and Uijk support this predietion. Examples are shown in Figures 1-3 for two salient 
non-equilibrium behaviours of the adaptive dynamies: a periodie attraetor and a ehaotie attraetor. 
In both eases, the figures show the eomparison of the adaptive dynamies ( [T3] ) with the eorrespond¬ 
ing individual-based model. Even though the individual-based simulation do of eourse not remain 
strietly monomorphie, they also do not break up into distinet elusters. Instead the phenotypie vari- 
anee remains eonstrained to a single eohesive eluster, and the trajeetory of the eenter of mass of 
this evolving eluster follows the adaptive dynamies trajeetory. Naturally, due to stoehastieity of 
the individual-based solution, the eenter of mass trajeetory is inevitably noisy. Nevertheless, the 
size, shape, and the direetion of motion of the individual-based trajeetory eorrespond to those of 
the adaptive dynamies solution. For ehaotie adaptive dynamies, this is a notable finding, as it was 
not at all elear a priori whether it is possible to find individual-based models that follow a ehaotie 
trajeetory without loss of eohesiveness and eoneomitant inerease of phenotypie varianee. 

A similar eomparison ean be performed for the adaptive dynamies and the eorresponding PDF 
model, although this is eurrently only feasible for dimensions d < 3. In sueh oases, the results 
are the same: the PDF model with the minimal eompetition kernel does not diversify, i.e., the 
phenotype distribution remains unimodal, and the movement of the mode of these distribution in 
phenotype spaoe traeks the trajeetory of the eorresponding adaptive models. This is illustrated in 
the right panels of Figures 1 and 2. 
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FIG. 1: An example of periodie adaptive dynamies in a 2-dimensional phenotype spaee (left 
panel), the trajeetory of the eenter of mass of the eorresponding minimal individual-based model 
(eenter panel and the first video (online only) in the seeond row ), and the trajeetory of the eenter 
of mass of the eorresponding minimal PDF model (right panel and the seeond video (online only) 
in the seeond row ) for a 2-dimensional system. Projeetions of the trajeetories onto the first two 
phenotypie eomponents are shown. The eoeffieients in the eompetition kernel ( [T?] ) are given in 

the file “eoeff2.dat” in SOM. 


B. Diversification with Gaussian competition kernels 


As mentioned before, the ehoiee of eompetition kernel Q for the given adaptive dynamies Q 
is by far not unique: Speeifieally, eonsider the set of eompetition kernels of the form 

d 


a(x,y) =exp 


^ “1“ Bijj^UiXjXk ~\~ Ciji^yiyjX]^ Diji^yiyjy]^ E^jXiXj H- FijyiXj -\- Gijyiyj) 


( 15 ) 
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FIG. 2: An example of ehaotie adaptive dynamies in a 3-dimensional phenotype spaee (left 
panel), the trajeetory of the eenter of mass of the eorresponding minimal individual-based model 
(eenter panel and the first video (online only) in the seeond row ), and the trajeetory of the eenter 
of mass of the eorresponding minimal PDF model (right panel, and the seeond video (online 
only) in the seeond row, only a small part of the attraetor is reprodueed sinee the PDE integration 
is very eomputationally extensive). Projeetions of the trajeetories onto the first two phenotypie 
eomponents are shown. The eoeffieients in the eompetition kernel ( [T4| ) are given in the file 

“eoef0.dat” in SOM. 


where the sets of eoeffieients {Aink}, {Bink}, {Cink}, {Anfc}, {Aj}, {Aj}, and {Gij} are subjeet 
to the following eonstraints: 


Aijk T Bijk -f Cijk Bijk 0 for all 

Eij + Fij -I- Gij = 0 for all i, j 


and 


(^ijk i^B^jk (^G{jk Gjik) T {^Dijk ~\~ Djik ~\~ Djki)^ 

^ij = “(Aj + {Gij + Gji)), 
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FIG. 3: An example of chaotic adaptive dynamics in a 4-dimensional phenotype space (left 
panel), the trajectory of the center of mass of the corresponding minimal individual-based model 
(right panel and the video (online only) in the second row ). Projections of the trajectories onto 
the first two phenotypic components are shown. The coefficients in the competition kernel ( [14] ) 

are given in the file “coeff4.dat” in SOM. 


where and bij are the coefficients determining the adaptive dynamics ( [T3j ). Then the first two 
of the above constraints ensure that a(x, x) = 1, and the third set of constrains ensures that the 
competition kernel ( [TS] ), together with the carrying capacity (j^, generates the adaptive dynamics 

0 . 

It is clear from these considerations that there are very many different competition kernels that 
give rise to the same adaptive dynamics. (In fact, the dimension of the space of order-3 competition 
kernels giving rise to one and the same adaptive dynamics model is 2cP + (f.) 

Here we are interested in those competition kernels corresponding to a given adaptive dynam¬ 
ics that would give rise to adaptive diversification. It has long been suggested that a general 
mechanism to generate and maintain diversity is for competition to be strongest between similar 


13 





phenotypes [|71 [TOl |24l|. This is typieally deseribed by Gaussian competition kernels, which re¬ 
flect the biologically plausible assumption that the strength of competition between individuals 
decreases with with phenotypic distance. 

In the present context, and starting with a given adaptive dynamics ( [T3| ), we can incorporate 
Gaussian competition by multiplying the minimal competition kernel Q with a Gaussian term, 
resulting in a competition kernel of the form 

/ W N [Xi — UiY 

^Wi[yi)[xi - Vi) - 
i=l 


a(x,y) = exp 


2a2 


(16) 


Note that this competition kernel is a particular case of the general form ( [Ts] ) (with {D} = 0 and 
a suitable choice of coefficients {A}, {B}, {C}, {E}, {F}, and {G}). Again, together with the 
carrying capacity ([^, this competition kernel results in the general adaptive dynamics ( fT3| ). But 
now the Hessian of the invasion fitness function has an additional positive diagonal component 


given by the cr. 


-2 


H(x) = - 


^ Wi — 


\Wd -Udj 


Wi - Ml, 


Wd 


Md) + 
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0 




-2 


a:(x) 
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This positive definite diagonal component is independent of the current resident phenotype, 
and hence in any phenotype space with dimension d > 2, there are directions in phenotype space 
orthogonal to the current selection gradient along which there is a disruptive component of se¬ 
lection. For sufficiently small variances ak of the Gaussian components, there will therefore be a 
tendency towards diversification. For any given point on an evolutionary trajectory, i.e., for any 
given resident phenotype, the conditions for diversification are of course not only determined by 


the positive definite component of the Hessian, but also by the remaining terms in (17) and the 
selection gradient (|7]). However, if the ak of the Gaussian components are sufficiently small, the 
diversification conditions are likely to persist for a sufficient amount of time to make diversification 
possible along non-equilibrium evolutionary trajectories. 

Whether evolutionary branching can occur under these conditions can be checked by using the 
adaptive dynamics ( [T3| ) from IfT^ as a starting point and constructing the corresponding individual- 
based model using the the competition kernel (16) and the carrying capacity ([^. For simplicity, 
we assumed that the variance of the Gaussian component ak = a was the same in all phenotypic 
directions, and our extensive simulations indicate that adaptive diversification indeed occurs when¬ 
ever a is small enough. Figure 4 illustrates how the typical interparticle separation, quantified as 


14 













FIG. 4: Time dependence of the total standard deviation of the particles’ coordinates defined as 
atotai = ~ 3-dimensional chaotic system as in Fig. 2 (bottom line), 

and with a Gaussian competition term ( [T^ with the width a = 0.65 (top line, red online). Here 
(...) stands for an average over all particles present in the system. Three snapshots for the 
Gaussian competition kerne illustrate the initial, single-cluster stage att = 40, intermediate 
diversification att = 1000, and the well-developed steady-state diversification att = 20000 are 
shown in inserts. The time evolution that has led to these snapshots can be seen in video in Fig. 6 


the “total standard deviation” atotai = y “ {xi)y), becomes noticeably larger for the 

individual-based model with the Gaussian competition kernel ( [T^ compared to the same model 
without the Gaussian term. The diversification in the models defined in Figs. 1-3 but with the 
Gaussian competition kernel is further illustrated by videos in Figs. 5-7. The corresponding PDF 
solutions also exhibit diversification, which is shown in videos in the right panels in Figs. 5,6. 

As mentioned, the Gaussian competition kernel ([T^ is a special case of the general competition 
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FIG. 5: (Online video) Diversifieation of the individual-based model (left panel) and the PDF 
solution (right panel) of the 2-dimensional periodie system presented in Fig. 1 but with a 


Gaussian eompetition term (161 with the width a = 0.65. 



FIG. 6: (Online video) Diversifieation of the individual-based model (left panel) and the PDE 
solution (right panel) of the 3-dimensional ehaotie system presented in Fig. 2 but with a Gaussian 


eompetition term (161 with the width a = 0.65. 



FIG. 7: (Online video) Diversifieation of the individual-based model of the 4-dimensional 


system presented in Fig. 3 but with a Gaussian eompetition term (16) with the width a 


ehaotie 
= 0.5. 
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kernel ( [T5| ), and there are potentially many more sets of eoeffieients in ( [T5| ) that generate the same 
adaptive dynamies ( [T3] ) but lead to diversifieation in the eorresponding individual-based models. 
Our goal here was not to enumerate all forms of a that ean produee evolutionary branehing, but 
rather to show that a partieular form of eompetition kernel generates models that do not diversify 
and follow the adaptive dynamies trajeetories, while eompetition kernels ineorporating suffieiently 
strong frequeney-dependenee leads to robust and persistent diversifieation in systems with non¬ 
stationary adaptive dynamies. 

IV. CONCLUSIONS 

Reeent advanees in understanding the evolution of eomplex multidimensional phenotypes were 
mainly limited to the adaptive dynamies framework IIT^ [131 |20l |32l or to the the quantitative ge- 
neties framework (blilSll. Neither of these frameworks lends itself direetly to the study of the 
evolutionary dynamies of diversifieation and subsequent eo-evolution of eoexisting phenotypie 
elusters, as adaptive dynamies assumes monomorphie populations, while quantitative geneties 
models typieally assume Gaussian phenotypie distributions. Here we have made an attempt to 
go beyond these approximations by studying stoehastie, individual-based model as well as deter- 
ministie PDE models for adaptive diversifieation under non-equilibrium evolutionary dynamies. 
We eonsidered a general elass of models for evolution due to frequeney-dependent eompetition. 
Starting with a given adaptive dynamies model in d-dimensional phenotype spaee (d > 1) with 
potentially eomplieated dynamies, we eonstrueted individual-based models that give rise to this 
adaptive dynamies model in the limit of large populations, and rare and small mutations. Without 
the assumption of small and rare mutations, the individual-based models give rise to determinis- 
tie PDE models in the large population limit ijSlSl, whieh in turn give rise to the given adaptive 
dynamies model in the limit of small varianee of the evolving phenotypie distributions. 

Our analysis generated two basie eonelusions: 

• Eor a given adaptive dynamies in multidimensional phenotype spaee there is a eorresponding 
“minimal” individual-based model whose fitness Hessian is always negative definite; as 
a eonsequenee, this minimal individual-based model does not diversify, and instead the 
population remains eonfined to a single eluster whose eentre of mass follows the trajeetory 
of the eorresponding adaptive dynamies model, regardless of the nature of the attraetor of 
the adaptive dynamies. 
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• Multiplying the minimal competition kernel by a Gaussian term with a sufficiently small 
width yields the same adaptive dynamics, but causes adaptive diversification in the 
individual-based model, in which multiple coexisting clusters emerge. Again this holds re¬ 
gardless of the nature of the attractor of the given adaptive dynamics. In particular, adaptive 
diversification is possible from a complicated evolutionary trajectory. 

Similar statements hold for the corresponding PDE models. In particular, diversification is 
possible in the absence of evolutionary equilibrium attractors in phenotype space, and in high¬ 
dimensional phenotype spaces diversification can in principle occur in any direction that is or¬ 
thogonal to the selection gradient. These findings considerably widen the scope of adaptive di¬ 
versification as a general evolutionary principle. We have also shown that the correspondence 
between individual;-based and PDE models on the one hand, and adaptive dynamics models on 
the other hand, is not unique: there are very many different individual-based and PDE models that 
give rise to the same canonical equation flSl for the monomorphic adaptive dynamics. Accord¬ 
ing to the analysis in [|2TI . diversification along complicated trajectories in these models should 
depend on the components of the Hessian matrix of second derivatives of the invasion fitness func¬ 
tion that are orthogonal to the selection gradients, and our results corroborate this. In particular, 
different “full” models reconstructed from a given adaptive dynamics have different orthogonal 
Hessians and different diversification properties. At one end of this set of full models are the min¬ 
imal models, which generate a negative definite orthogonal Hessian, and therefore do not show 
diversification. The other extreme is the Gaussian model, which can be defined to have Hessians 
with positive eigenvalues in all orthogonal directions, regardless of the current resident phenotype. 
Consequently, these reconstructed “Gaussian” models have a high propensity to diversify. 

It is interesting to note that the minimal models reproduce the adaptive dynamics attractor even 
if this attractor is chaotic. Due to the intrinsic sensitivity and complexity of chaotic systems, it 
was not clear a priori that this is possible. Eor the Gaussian models presented here we made 
the simplifying assumption that competition between similar phenotypes is equally strong in all 
phenotypic directions. In reality, Gaussian competition could act only in a subset of phenotypic 
directions, and then one would expect diversification to primarily occur in this subset. Also, we 
assumed all off-diagonal elements in the Hessian portion of the competition kernel to be 0, which 
corresponds to assuming that there are no interactions between the phenotypic components with 
regard to the Gaussian component of competition. It is known that such interactions promote 
diversification in equilibrium models [|6l[I3|32l, and it is an interesting direction for future work 
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to analyze whether similar results hold for diversification in non-equilibrium models. 

Our simulation results indicate that diversification often leads to coexistence of multiple pheno¬ 
typic clusters, and a very interesting question for future work concerns the effect of diversification 
on the complexity of the co-evolutionary dynamics of these coexisting clusters. So far, we have 
seen examples of both destabilization and stabilization due to adaptive diversification: popula¬ 
tions converging to an equilibrium when monomorphic can embark on non-equilibrium dynamics 
after diversification; and the opposite can happen as well, so that populations moving on a com¬ 
plicated attractor when monomorphic converge to a multi-cluster evolutionary equilibrium after 
diversification. 

It seems clear that in general, many different phenotypic properties can contribute to ecological 
interactions, resulting in potentially complicated selection pressures in high-dimensional pheno¬ 
type spaces. Our previous work has highlighted both the increased propensity for evolutionary 
branching in such spaces [IT3]I . and the increased propensity towards complicated evolutionary 
trajectories IfT^ . Here we have extended this perspective to studying the full dynamics of diver¬ 
sification in high-dimensional phenotype spaces in stochastic, individual-based models as well as 
PDE models. When these evolving systems diversify, the resulting coevolutionary dynamics of co¬ 
existing phenotypic clusters may often be even more complicated and unpredictable, thus further 
challenging the prevailing view of evolution as an optimization and equilibrium process. We hope 
that our work will contribute to the discussion about the importance of adaptive diversification as 
a mechanism for generating biological diversity. Finally, we think that studying individual-based 
and PDE models of diversification in high dimensional spaces can also be relevant for tackling 
important questions in cultural evolution, such as the origin and evolution of different religions, 
languages, and other cultural traditions lITOl . 
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Appendix A: Supplementary online materials 


Here we deseribe how we simulate the individual-based and PDE models introdueed in the 
main text. Both models are defined by the eompetition kernel 



d d / \2 


(Al) 


and the earrying eapaeity 



(A2) 


Sinee we are mainly interested in non-equilibrium dynamies, the sets of eoeffieients {a} and {6} 
and the initial eonditions xq were seleeted as deseribed in llT^ . so that the eorresponding adaptive 
dynamies, given by Eq. 13 in the Main Text, was either eyelie, quasiperiodie or ehaotie. 

1. Individual-based simulation 

Individual-based realizations of the model were based on the Gillespie algorithm ifTTll and eon- 
sisted of the following steps: 

1. The system is initialized by ereating a set of Kq ~ 10^ — 10'^ individuals with phenotypes 
Xq, G R'^ loealized around the initial position xq with a small random spread | Xq — xq | ~ 


2. Eaeh individual a has a eonstant reproduetion rate pa = 1 and a death rate 6a = 

A(xq, X ij)/K{yia), as defined by logistie eeologieal dynamies. 

3. The total update rate is given by the sum of all individual rates U = + ^a)- 

4. The running time t is ineremented by a random number At drawn from the exponential 
distribution P{At) = U exp{—AtU). 

5. A partieular birth or death event is randomly ehosen with probability equal to the rate of this 
event divided by the total update rate U. If a reproduetion event is ehosen, the phenotype 
of an offspring is offset from the aneestral one by a small mutation randomly drawn from a 
uniform distribution with amplitude /i = 10“^ — 10“^. 
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6. The individual death rates da and the total update rate U are updated to take into aeeount the 
addition or removal of an individual. 

7. Steps 4-6 are repeated until t reaehes a speeified end time. 


2. Partial differential equation models 


A deterministic large-population limit of the individual-based model is obtained as the partial 
differential equation (PDE) 


aiV(x, t) ^ ^ 1 - f a(y,x)N(y,t)dy \ ^ ^ ^ d^N(x, t) 


dt 


K{x) 


i=l 


dxj 


(A3) 


where N{x,t) is the population distribution at time t [|3l. The second term of the right hand 
side is a diffusion term that describes mutations, with the diffusion coefficient typically set to 
D ~ 10“^ — 10“^. The form and size of the single-cluster trajectory was usually known from 
the adaptive dynamics solution. Hence, to numerically solve the PDE model ( |A3| ) we chose a 
finite lattice resolution of phenotype space nodes at least twice larger than the adaptive dynamics 
attractor. The number of bins L in each dimension of this lattice is strongly limited by memory 
limitations: An efficient implementation requires computing and storing an array of x values 
of the competition kernel Q;(yj, xj) for the pairwise interactions between all sites i and j. With 
L = 25 — 30 to achieve a reasonable spatial resolution, the memory constraint makes the PDE 
implementation feasible only for d = 2,3. 
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